arXiv: 1501.02168vl [math.NA] 9 Jan 2015 


Convergence and visualization of 
Laguerre’s rootfinding algorithm 

Herbert Moller* 


Version of January 9, 2015 

Abstract Laguerre’s rootfinding algorithm is highly recommended although 
most of its properties are known only by empirical evidence. In view of this, we 
prove the first sufficient convergence criterion. It is applicable to simple roots 
of polynomials with degree greater than 3. The “Sums of Powers Algorithm” 
(SPA), which is a reliable iterative rootfinding method, can be used to fulfill the 
condition for each root. Therefore, Laguerre’s method together with the SPA is 
now a reliable algorithm (LaSPA). In computational mathematics these results 
solve a central task which was first attacked by L. Enler 266 years ago. 

In order to study convergence properties, we eliminate the polynomial and its 
derivatives in the definition of the Laguerre iteration, replacing them by sums, 
only depending on the roots and the iterated values. For this iteration with roots, 
the above criterion of convergence represents an efficient stopping condition. In 
this way, we visualize convergence properties by coloring small neighbourhoods 
of each starting point in squares. 
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1. The convergence criterion 

If q is a simple root of the polynomial p(z), then it has been proved that Laguerre’s 
method converges cubically to the limit g, whenever the initial guess is close 
enough to g. But it was not known how small the distance from the root must 
be. In the following theorem we determine for each simple root g a disk with 
centre g such that for all starting values in this neighbourhood, the convergence 
to the root is guaranteed. If all roots of p(z) are simple, an a priori lower bound 
for the radii of all these disks is given. At the end of this section, we explain how 
the condition of the theorem can be reliably fulfilled with the SPA. 

Theorem. Let p(z) be a normalized polynomial with degree m > 4 and let P : = 
{(Oi, ..., Qi} be the set of the roots. For w 0 e C \ P with p'(uo) ^ 0 or p"{uq) ^ 0 
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the Laguerre sequence L p (u 0 ) =: ( u n ) n£ ^ is defined recursively by 


m 


Mn+l — U r 


n G N, with q(z ) : = 


p'(z) 

q{u n ) + s(u n )r(u n y —’ 

P " W and 


(1) r(z) := \J(m- l)(m((z) - j 2 (z)), «(z) : = q 2 (z) 


p(z) 


s z : = 


1 when (Reg(z)) (Rer(z)) + (lmg(z))(lmr(z »>°. 
-1 else, 


as long as q(u n ) + s(u n ) r(u n ) 0. 

If there is a simple root g G P such that 


( 2 ) 


K - Q\ < 


2 m — 1 


nun 


in {Io’ - Q\ | cr E P\ {ft}-}, 


then L p (u 0 ) converges to the limit g, and 


( 3 ) 


15 

\u n ~ q\ < A n \u 0 — with A : = — 

16 


holds for all n G N \ {0}. 

If all roots of p(z) are simple, then in (2) p, Q : = min {|cr — £>| | cr G P \M} can 
be replaced by the a priori lower bound 


1 + ■p-r max ^ \dj\) 2 < minj/ig | g G P}, 


i<j<( 


where dj are the coefficients of the polynomial D p (z) : = (^ — (ft — Pk ) 2 ) = : 

l<z<fc<m 

(?) 

Y djzf These coefficients can be calculated only using the coefficients of p(z). 
l=o 

All roots of p(z) are simple if and only if d 0 0. 


Proof. From p(z) = n i z ~ Qk) mk , where rn k is the multiplicity of g k for k = 
k =1 

1,... ,1, it follows that 


q{z) = ZJZV = w-1 np(z) = ^^m fc ln (z - p k ) = 


( 4 ) 


t(z) = q\z) 


p{z) dz 
p"(z) 


k =1 


k =1 


m k 

z- Pk 


P(z) 


dz 




m k 


Y ( z ~ ^’) 2 


and 
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Without loss of generality, g can be chosen as Q\ with m i = 1. With the abbre¬ 
viations 


( 5 ) 


Si = (n 0 ) : = —, S 2 = S 2 (u 0 ) : = V 

^' 7/o — /It ^' 


k=2 


Uq ~ Qk 


= 2_^ 

k=2 


Uq ~ Qk 


and 


w = w(uo) : = ^ (l + Si + s \J{m — l)(m(l + S 2 ) ~ (1 + Si) 2 ) ), 
where s = s(n 0 ) G {—1,1} is defined in the same way as s(u 0 ), we get 

l«i - el = l«o - el I 1 - R 


( 6 ) 

Since the condition (Req(z)) (Rer(^)) + (lmg(z))(lmr(z ))> 0 is equivalent to 
| q(z) + r(z) | > | q(z) — r(z) | and since the latter relation doesn’t change when 
each term is multiplied by the same non-zero factor, it follows that s(uo) = s(uo)- 
First, it will be shown that (2) implies s(w 0 ) = 1. 

From (2) we have (2m — l)|it 0 — £?| < W~ q\ — |cr — uq + Uq — < \uq — a| + |no — q\ 

for each a G P \ {£>}, from which I——-I < „ 1 „ follows. Therefore and with 


u 0 -a i 


2 m —2 


^2 rrik = m — 1 we get 

k=2 

( 7 ) 


N < k 1*1 < id 


4m—4' 


With Sj : = Re Sj and tj : = Im Sj, j = 1, 2, (7) leads to 


( 8 ) 


|sr| < L \ti\ < l, |s 2 | < 


1 


4m—4 


and \t 2 \ < 


1 


4m—4' 


( 9 ) 


We evaluate the square root in (5) with the well-known formula 

Va + ib = + a) + i (sign b)b 2 - a) 

for a, b G M and 6^0. 

Since r : = \J(m — 1)(m(1 + S 2 ) — (1 + Si) 2 ) = : y/(m— l)(c + id) with c = 
m — 1 + ms 2 — 2Si — + t\ and d = mt 2 — 2 ti — 2sRi, (9) gives 

Re r = \J -^ 2 ^- (V c 2 + d 2 + c) and 

Im r = (sign d) \J (Vc 2 + d 2 — c). 

qi n 11 

From 0<m—j^<c<m+j^ and \d\ < for m > 4, we get 

yj(m — 1 )(m — < Re r < yJ (m — 1) (m + |) < m + ^ and 

(10) 


Im rl = / 11 

2(V^TdP + c) - 12' 


m — 1 11 

m - v/68 
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1 r* A Q 1 O 7 

With the abbreviations a g : = ^ := ; starting from 

?n > (7 + \Jg 2 — h = 3.99639 ... and passing m 2 — 2gm > — h, we get from (10) 

( 11 ) 


16 1 

Rer > —m -for m > 4. 

25 2 


The estimates ( 8 ), (10) and (11) lead to (Re (1+S'i)) (Re r) + (lm (1+Si)) (im r) > 


2(25™ ~ 2) 


2-\/68 


> 1 for m > 4, which yields 
s(m 0 ) = s(u 0 ) = 1 . 


Next, continuing with ( 6 ), we calculate an upper bound for d : = |l 
a; : = Re«j and y : = Inno. Since 

(12) * = =,/Ti-^Wr + 




with 


a ; 2 + a / 2 


a : 2 + a / 2 


2 / 


x 2 + a / 2 


2 = ,/i- 21 


x 2 + y- 


we need a lower bound for x and an upper bound for x 2 +y 2 . From ( 8 ), (10) and 
( 11 ), we obtain 

x = R 1 + + Ref ) > R 1 -1 + i m - 5 ) = i’ 


m K " ' m ' 2 

+ \ + m + ^ i and 


M < ~(ti + |Im r|) < ^- + _) < - for m > 4. 


Therefore with (12), we get 


(13) 


2162577 ^ 15 _ 
2465425 < 16 ~~ 


Now, we prove (3) by mathematical induction. With n — 0 and 11 — ~ 
( 6 ) and (13) result in the basis of the induction 


a? < A, 


K - p| < A |w 0 - g\. 

As inductive step, we show that \uh — q\ < X h |u 0 — £?| leads to \uh+i — g\ < 
X h+1 |w 0 — p|, if h is a generic number. Since X h < 1, using ( 2 ), we have 

(14) \u h - g\ < X h \u 0 - g\ < 1 \x e . 

Therefore, in the above derivation, we may replace u 0 by Uh and u\ by Uh+ 1 - Then 
we get 

\uh+i ~ g\ < X \u h - g\ < X h+1 | u 0 - Pi- 

Reasoning by induction, it follows that (3) is valid for all n G N \ {0}, which 
means that L p {uq ) converges to the limit g. 
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The a priori lower bound for the distances of the roots of p(z) is derived in [3] (p. 
379). The coefficients of D p (z ) are calculated with the aid of the sums of powers 

m 

°j ■ = 3 = !, • • •, ( m - l ) m and <r'j : = ( Qi - 3 = !, ■ ■ ■, (• 

k= 1 1 <i<k<m 

771—1 

With the coefficients of p(z) = : Y CkZ k + z m , we have 

k =0 

{ k—1 

— Y Cm-j&k-j ~ kc m _k for 1 < k < m and 

m 1 

— Y2 c m-j&k-j for m < k < (m — 1 )m. 

j =i 

Expanding the powers in crj, we get 

o' = ma 2j + J2(- 1 ) l ( 2 l) ai a V- 1 + ( _1 ) J { 2 j-i) a j for j = 1, • • •, (™) • 

i=i 

Reversing the first case of the above recursion for cr^ with a' k instead of a k and 
replacing c m _j with dtm\_j, we obtain 

k—l 

d (™)-k = d (™)-j a k-3 fOT k = • • • > (T) • 

j=0 

The lower one of the estimates in [3] (p. 375, formula (100)) applied to D p (z ) 
gives 

(1 + max IcLI ) 2 < min{u 0 I p G Pj. 

V l d °l t<i<(™) / e 

Since D p ( 0) = 0 if and only if d 0 = 0, all roots of p(z) are different if d 0 ^ 0. □ 

These sufficient criteria for the convergence of Laguerre sequences can be reliably 
fulhlled for each root with the aid of the SPA, which was introduced in the book 
[3]. A revised version was published in [4] and [5]. Therefore we will only sketch 
how the reliability is gained and which of the used methods are also valuable for 
Laguerre’s method. 

The SPA starts at a present version of D. Bernoulli’s method which uses quotients 
of successive sums of powers of the roots to approximate dominant roots. We 
have proved two new basic results for the sequence of the quotients which we 
call “Bernoulli sequence”, namely, a recursion formula and a representation of 
the minimal absolute value of the roots only using the terms of the Bernoulli 
sequence. 
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In the case of (nearly) equimodular roots, the second result is combined with the 
search for local minima of the absolute values of the polynomial with arguments 
on a circle which has a good approximation of the minimal modulus as radius. 
All these minima are close to roots, and at least one of these roots lies inside the 
circle which we call “minimum circle”. 

If necessary, this construction is repeated with “chained minimum circles” or with 
“modified Turan circles” which have radii forming a null sequence. Therefore, in 
any case, after finitely many steps a root can be approximated with prescribed 
accuracy. 

Error bounds and stopping criteria are obtained with “Laguerre disks” [1] 

C u : = {w € C | \w — u + 7Tip(u)/(2p'(u)) [ < m \ p(u) / (2 p'(u)) |} 
for u £ C with p(u) 7 ^ 0. 

Each Laguerre disk contains at least one root g £ P, the radius of £„ is smaller 
than m\u — £>| if 2m\u — £>| < min{/J e | g £ P}, and in this case, there is only one 
root in C u . With Laguerre disks the SPA reliably separates all roots in disjoint 
disks. Using the distances of such disks as lower bounds for p g in (2), the above 
Theorem guarantees that the SPA can be terminated with Laguerre’s method. 

Up to now, computer programs for Laguerre’s method use values of the polyno¬ 
mial in the stopping criteria (see, for example, the NAG routine presented in [ 2 ] 
and the C program in [9]), and, if at all, there are only complicated error bounds. 
Now, with the above Theorem and with Laguerre disks, we have precise stopping 
criteria and error bounds also for Laguerre’s method. 

If a preset bound for the step numbers is exceeded using the very good but un¬ 
proved convergence properties visualized in the next section, then the SPA will 
serve as a “security net”. In this way, Laguerre’s method together with the SPA 
is not only highly efficient but also reliable. Therefore we call this combination 
“Laguerre and SP Algorithm” (LaSPA). The connection with Euler’s investiga¬ 
tions in his famous book “Introductio in Analysin Inhnitorum” (Introduction to 
Analysis of the Infinite) is explained in [4] and [5]. 

2. Visualization of convergence properties 

The fundamental theorem of algebra constitutes a one-to-one correspondance 
between finite sets of complex numbers and normalized polynomials with simple 
roots. Therefore, in (1) we may use the representations of q(z ) and t(z ) by the 
terminating sums in (4) to get Laguerre sequences only depending on the roots 
and the iterated values. Since, in that way, it is very easy to fulfill the condition of 
(2), we have written a Cython program “Laguerre.pyx” to visualize convergence 
properties of Laguerre’s method in squares which can be freely chosen parallel 
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to the axes in the complex plane. Onr first goal was to check the statement of 
the following sentence in the introductory paragraph of the lemma “Laguerre’s 
method” in the English part of Wikipedia: <One of the most useful properties 
of this method is that it is, from extensive empirical study, very close to being a 
“sure-fire” method, meaning that it is almost guaranteed to always converge to 
some root of the polynomial, no matter what initial guess is chosen.> 

The program consists of 150 lines. We will only shortly explain it, because it 
is available in the section “English subjects” of [8], and it is commented in [7]. 
The Sage system [10] is needed to run the program, because Laguerre.pyx uses 
graphics from Sage, and Cython, after being translated to C, is compiled by the 
GNU C compiler included in Sage. 

The program takes six items as input: a list of different complex numbers (the 
roots), the midpoint of a non-rotated square in the complex plane, the length of 
the sides of this square, the number m of starting values in each row and column 
with equal distances inside the square, a number for the image to be saved, and 
the extension of the image file (eps, pdf, png and 3 more). If 0 is entered as the 
length of the sides, a “standard square” will be used which has the same midpoint 
as the smallest non-rotated rectangle containing all roots, and the length of the 
sides is set to the double of the maximal length of the sides of the rectangle. 

The program mainly consists of three parts. At the beginning, most of the 
parameters get their initial values, and the squared minimal distance of each 
root from the other roots is calculated. In the central part, for each of the m 2 
starting values the terms of the corresponding Laguerre sequence are computed 
until one of the conditions in (2) is fulfilled or a preset bound for the step number 
is succeeded. In the first case, the starting value is stored in two lists, namely 
one belonging to the root determined by (2), and the other assigned to the step 
number. In the second case, the starting value is recorded as an indication for a 
cyclic Laguerre sequence. 

The last part organizes the output which consists of two figures and a list, report¬ 
ing for each occurring step number how many starting values need this number 
of steps to fulfill the condition of (2). The figures will be explained and analysed 
using an example with five roots and 250000 starting values which we will call 
“starting points” in the following. 

Both figures contain five white areas which represent the disks determined by (2). 
We call each of these disks “safe”, because a Laguerre sequence with a term in 
such an area converges safely to the root in the centre of the disk. The list of the 
roots is [1.6 - 0.55.7, -0.39 + 0.03j, -2.32 + 2.17j, 0.2 - 1.06j, -0.02 - 0.27j], 
where the terminating j or J is the notation of Cython for the imaginary part 
of complex numbers. All starting points in the rest of the figures have coloured 
neighbourhoods which are touching disks in the case of PDF files and small 
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squares with shading colours for the pixel graphics of PNG output which can be 
converted to EPS format. 

On a computer with a 2.4 GHz processor, the time for translating to C and 
compiling was 10 seconds, the data were calculated in 8 seconds, and after 110 
seconds the graphics for two times the starting points outside the safes appeared. 
Since with 90000 starting points the latter time was 50 seconds, the program can 
also be used to produce short motion picture sequences, for example, to visualize 
the effect of moving zeros or for zooming. 

The coloured areas in Figure 1 are named “limit areas”, because the neighbour¬ 
hood of the starting point of a convergent Laguerre sequence gets the colour of 
the area surrounding the safe which contains the limit of the sequence. It is one of 
the hypotheses formulated afterwards, that these points indeed are not isolated. 



Figure 1. Safes (white disks) and limit areas 
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In FIGURE 2, the neighbourhood of the same starting points as in FIGURE 1 
is coloured according to the number of steps until the corresponding Laguerre 
sequence reaches a safe. Therefore, the areas with the same colour are called 
“step areas”. If the occurring colours are listed in the order of the rainbow from 
red to violet, then the index of the colour is the step number. For example, the 
area surrounding a safe is always red, because one step is needed to enter the 
safe. Moreover, the second figure shows a trajectory from a starting point to a 
safe with the maximal number of steps, if the trajectory completely lies in the 
square. 

The arithmetic mean of the roots is given by a : = ——c m _i. If it fulfills V -^ 

& j m m 1 ^ a-Qk ' 

m m 

0 or V - + V t- 77 , which means that q(a) + s(a) rfa) + 0 , and if it lies 

^ a—Qk / ^ {a-BkY \ \ > I ’ 



Figure 2. Safes, step areas, a trajectory with maximal step num¬ 
ber 7, and the arithmetic mean of the roots (black point) 
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in the square, then it is shown in the second figure as a black point, because in 
most cases it is an optimal starting point. 

The possibility of selecting arbitrary squares turns the program into a tool. But 
here, having compared many different figures, we can only state several hypothe¬ 
ses and a conjecture concerning the convergence properties of Laguerre’s method. 
The corresponding figures with comments and the description of modified pro¬ 
grams are contained in a report [7]. 

• The complex numbers as possible starting points can be divided into four 
types: inner points, boundary points, “cycle points” which generate cyclic 
Laguerre sequences and “singularity points” which as starting points are 
excluded by definition, because they are the multiple zeros of p'(z). 

• The inner points and the boundary points constitute the limit areas and 
the step areas which both have boundaries consisting of piecewise smooth 
curves. If there are no singularity points, then a number S p exists such 
that the starting points of same type with modulus greater than S p form 
unbounded and simply connected areas. Moreover, each corresponding 
part of the boundary curve of a limit area or a step area is asymptotic to 
a straight line. 

• The cycle points and the singularity points are the only isolated points. 

• The boundary curve of each step area belongs completely or piecewise to 
one of the neighbouring areas. 

• If C is a singularity point or a cycle point, then for each number Af 6 N 
there exists a disk A ^ with centre ( such that the step numbers of all 
Laguerre sequences with u 0 G A N \ {£} are greater than N. 

• All Laguerre sequences are bounded. This shall be a conjecture, because, 
possibly, it can be proved with a method similar to that used for the 
Theorem above. Namely, in this way, it is easy to show that A tends to 
0 if Uo increases unboundedly. 

One of the modified programs verifies that it is important to use the sign in 
the denominator of Laguerre’s method, because without this condition the limit 
areas and the step areas become highly fragmented. A second program is de¬ 
signed for the positioning of very small squares using the list with the occur¬ 
ring step numbers and turning off the graphics. In this way, for the polynomial 

p(z) = z 5 + (0.93 - 0.32 i)z 4 - (0.5818 - 5.8351 i)z 3 - (6.250301 + 1.40811 i)z 2 - 
(1.61432911 + 2.57707253 i)z + 0.306917325 - 0.5317505851, which corresponds 
to the above root list, the (hypothetical) cycle (0.063005... + 0.1196569... i, - 
0.327219... - 0.4305399... i) of length 2 is determined. 
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